#Load required libraries

library(tidyverse)
library(gsw)
library(ggplot2)
library(plotly)
#install.packages('seacarb')
#install.packages('prettydoc')

#Now we need to import our data

library(readr)
hydrostation_bottle <- read_table("hydrostation_bottle.txt", 
    col_names = FALSE, skip = 31)

library(readr)
hydrostation_bottle_names <- read_csv("hydrostation_bottle.txt", 
    skip = 30)


colnames(hydrostation_bottle)=colnames(hydrostation_bottle_names)
#view(hydrostation_bottle)

Variable Names and Units

#lets firts plot the data
hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x = decy, y = `Sig-th`))

hydrostation_bottle %>%
  filter(`Sig-th`!=-999) %>%
  ggplot()+
  geom_point(aes(x = decy, y = `Sig-th`))

hydrostation_bottle %>%
  filter(`Sig-th`!=-999 & Depth <20) %>%
  ggplot()+
  geom_point(aes(x = decy, y = `Sig-th`))

hydrostation_bottle %>%
  filter(`Sig-th`!=-999 & Depth <20) %>%
  ggplot()+
  geom_line(aes(x = decy, y = `Sig-th`))

#clear seasonal signal for sigma-theta,lets see how this compares to temperature
hydrostation_bottle %>%
  filter(`Sig-th`!=-999 & Depth <20) %>%
  ggplot()+
  geom_point(aes(x = Temp, y = `Sig-th`))

#Temperatures and density ares strongly correleted, but there appears to be two outliers tht we sill likley need to adreess at some point
#we only have density data from 1988-present, but temperature and salinity data from 1950s-present 
#This means I can calculate seawater density from 1950s-present 
#we can use the TEOS-10 to do this

TEOS-10 Toolbox in Package seacarb

?gsw #launches the documentation for the gibbs seawater toolbox (TEOS-10). All of this function within the package 

?gsw_sigma0 #lets check this function 
#it says we need absolute salinity and conservative temperature USAGE

#first we need absolute salinity:
?gsw_SA_from_SP
#practical salinity
#sea pressure (dbar)
#longitude 
#latitude

#Lets plot our pressure data - its missing before 1980s
hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=Pres))

#we have depth data for the time series 
hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=Depth))

?gsw_p_from_z

#adds a pressure colum from the depth and latN colums from/to hydrostation bottle #creates new variable

hydrostation_bottle=
hydrostation_bottle %>%
  mutate(Pres_gsw=gsw_p_from_z(Depth*-1,latN))

hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=Pres, y=Pres_gsw))
## Warning: Removed 3 rows containing missing values (`geom_point()`).

#we see strong 1:1 agreement between measured pressure and calculated pressure

hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=latN))

hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=lonW))

hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=CTD_S))

#Checking lat, lon and salinity data
hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=Sal1))

hydrostation_bottle=
hydrostation_bottle %>%
  mutate(Pres_gsw=gsw_p_from_z(Depth*-1,latN)) %>%
  mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN))

#chek it 
hydrostation_bottle %>%
  ggplot()+
  geom_point(aes(x=decy, y=S_abs_gsw))
## Warning: Removed 3 rows containing missing values (`geom_point()`).

#how else can I check my data
hydrostation_bottle %>%
  filter(Sal1!=-999) %>%
  ggplot()+
  geom_point(aes(x=Sal1, y=S_abs_gsw))

#now we need to calculate conservative temperature
?gws_CT_from_t
## No documentation for 'gws_CT_from_t' in specified packages and libraries:
## you could try '??gws_CT_from_t'
#we need absolute salinity, in-situ temp (ITS-90), and sea pressure 

#add line to calculate conservative temperature and change from hydrostation_bottle to HydrosS to not create confusion
HydroS=
hydrostation_bottle %>%
  filter(Sal1!=-999) %>% #filter for salinity missing
  filter(Temp!=-999) %>%
  mutate(Pres_gsw=gsw_p_from_z(Depth*-1,latN)) %>% #pressure for depth
  mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN)) %>%
  mutate(T_cons_gsw=gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw))

#missing data por eso se añadio filter 

#hydrostation_bottle %>%
  #filter(Sal1!=-999) %>%
  #ggplot()+
  #geom_point(aes(x=Sal1, y=S_abs_gsw))

#hydrostation_bottle %>%
  #filter(Temp!=-999) %>%
  #ggplot()+
  #geom_point(aes(x=Temp, y=T_cons_gsw))

HydroS=
hydrostation_bottle %>%
  filter(Sal1!=-999) %>% 
  filter(Temp!=-999) %>%
  mutate(Pres_gsw=gsw_p_from_z(Depth*-1,latN)) %>% 
  mutate(S_abs_gsw=gsw_SA_from_SP(Sal1,Pres_gsw,360-lonW,latN)) %>%
  mutate(T_cons_gsw=gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw)) %>%
  mutate(sig_th_gsw=gsw_sigma0(S_abs_gsw,T_cons_gsw))
  
HydroS %>%
  filter(`Sig-th` !=-999) %>%
  ggplot()+
  geom_point(aes(x=`Sig-th`,y=sig_th_gsw))

#but we have a very low sig-th-gsw so lets find it 

#our bottle and ctd salonity do not agree, this is likley the problem for the low
HydroS %>%
  filter(sig_th_gsw<0) %>%
  view()

HydroS_correctedS_a=
  HydroS %>% 
  filter(sig_th_gsw<0) %>% #filter for negative value 
  mutate(S_abs_gsw=gsw_SA_from_SP(CTD_S,Pres_gsw,360-lonW,latN)) %>%
  mutate(T_cons_gsw=gsw_CT_from_t(S_abs_gsw,Temp,Pres_gsw)) %>%
  mutate(sig_th_gsw=gsw_sigma0(S_abs_gsw,T_cons_gsw)) 

HydroS_correctedS_b=
  HydroS %>% 
  filter(sig_th_gsw>0) 

HydroS_correctedS = rbind(HydroS_correctedS_a,HydroS_correctedS_b)

HydroS_correctedS %>% 
  filter(`Sig-th`!=-999) %>%
  ggplot()+
  geom_point(aes(x=`Sig-th`,y=sig_th_gsw))

#Homework 

HydroS_correctedS %>% 
  ggplot()+
  geom_point(aes(x=sig_th_gsw,y=Depth))+
  scale_y_reverse()+
  scale_x_continuous(position='top')+
  #xlab(expression(paste(sigma[theta], "(kg m"^"-3",")"))) 
  xlab(expression(paste(sigma*theta, "(kg m"^"-3",")" )))+
  ylab('Depth (m)')+ 
  theme_classic()

Has surface sigma theta decreased over time?

HydroS_shallow= HydroS_correctedS %>%
  filter(Depth>30)

?lm #constroc a linear model
#lm(y~x, data=data) y is a funcion of x meaning data=data 
lm(sig_th_gsw~decy,data=HydroS_shallow)
## 
## Call:
## lm(formula = sig_th_gsw ~ decy, data = HydroS_shallow)
## 
## Coefficients:
## (Intercept)         decy  
##   2.587e+01    5.271e-04
#coefficients (intercept and decay)
#y=mx+b
#y=sig_th_gsw
#x=decy
#coefficient: intercept= b, decy=m
#sig_th_gsw= -0.004*decy + 33.4
#(kg/m3) = (kg/m3/y)*y + (kg/m3)

sig_theta_time_model=lm(sig_th_gsw~decy,data=HydroS_shallow)
summary(sig_theta_time_model)
## 
## Call:
## lm(formula = sig_th_gsw ~ decy, data = HydroS_shallow)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4145 -0.5393 -0.1192  0.8283  1.7717 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.587e+01  4.779e-01  54.130   <2e-16 ***
## decy        5.271e-04  2.402e-04   2.195   0.0282 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7605 on 26275 degrees of freedom
## Multiple R-squared:  0.0001833,  Adjusted R-squared:  0.0001452 
## F-statistic: 4.816 on 1 and 26275 DF,  p-value: 0.0282
plot=HydroS_shallow %>%
  ggplot(aes(x=decy,y=sig_th_gsw))+
  geom_point()+
  geom_line()+
  geom_smooth(method="lm")+
  theme_classic()

ggplotly(plot)
## `geom_smooth()` using formula = 'y ~ x'

#Lab Assignment

  1. Pick a question
  2. Produce a plot and a statistical summary using lm()
  3. Describe your results, the summary, and answer the question
  4. Compile into a completed lab report using R markdown

Potentia question: Include hypothesis (usar pag de cosas que se pueden calcualr)

How do temperature, salinity, and sigma-theta co-varyÇ Is there a relationship between sigma theta and oxygen? Is there a relationship of the parameters with depth? with time? Whithin a depth range over time? Are there seasonal differences in any of the parameters?